/* norms.h
 *
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#pragma warning(disable:4996)

#if !defined(_norms_)
#define _norms_
#include<vector>
#include <src/miscellaneous/OS_basic.h>
#include <src/data/Matrix.h>
#include <cmath>
/*! \file norm.h
	 \brief tools to define norms of multidimensional vectors 
*/
	  namespace jmmath {
			template<class T, class B>
			struct assign_helper {
				inline void operator()(T& t, typename T::base_value_type val)
				{
					t.assign(val);
				};
			};
			template<class T>
			struct assign_helper<T,OS_true> {
				inline void operator()(T& t, T val)
				{
					t = val;
				};
			};

/*! 	 \brief perform the sum \f$ \sum_i a_i \f$
*/
		  template <typename T >
		  typename T::value_type sum( const T & Array) 
		  {
			  typedef T::value_type value_type;
			  if (Array.size() == 0) return T::value_type();
			  value_type local(*Array.begin());
			  assign_helper<value_type, boost::is_fundamental<value_type>::type>()(local,0);
			  T::const_iterator it_b = Array.begin(),it_e = Array.end();
			  for (; it_b != it_e; ++it_b) {
				  local += (*it_b);
			  }
			  return local ;
		  };
/*! 	 \brief returns \f$ \sup_i a_i \f$
*/
		  template <typename T >
		  OS_double norm_infty( const T & Array) 
		  {
			  if (Array.size() == 0) return 0.;
			  OS_double local(0.);
			  local = 0.;
			  T::const_iterator it_b = Array.begin(),it_e = Array.end();
			  for (; it_b != it_e; ++it_b) {
				  local = std::max( norm_infty(*it_b), local) ;
			  }
			  return local ;
		  };
/*! 	 \brief returns \f$ sum(a*b)  \f$, that is usually a scalar product.
*/
		  template <typename T1, typename T2 >
		  OS_double scalar_prod( const T1 & p1, const T2 & p2) 
		  {
			  return sum(p1 * p2);
		  };
/*! 	 \brief returns \f$ \frac {sum(a*b)}{a.size()}  \f$, that is usually a normed scalar product.
*/
		  template <typename T1, typename T2 >
		  OS_double normed_scalar_prod( const T1 & p1, const T2 & p2) 
		  {
			  return sum(p1 * p2)/p1.size();
		  };

/*! 	 \brief returns \f$ \frac { \sum_i a_i }{a.size()}  \f$.
*/
		  template <typename T >
		  typename T::value_type moment1( const T & Array) 
		  {
			  return sum(Array) / Array.size();
		  };
/*! 	 \brief returns \f$ \frac { a*a }{a.size()}  \f$.
*/
		  template <typename T >
		  typename OS_double moment2( const T & Array) 
		  {
			  return sum(Array*Array) / ( Array.size() );
		  };
/*! 	\brief Variance calculus.
		Variance calculus for the equiprobable case : \f$ /frac{1}{N}\sum_{i=0,..,N-1} |x_i|^2 - | \overline{x} |^2 \f$
*/
		  template <typename T >
		  typename OS_double Var( const T & Array) 
		  {
			  return moment2(Array) - norm22( moment1( Array) );
		  };
/*! 	\brief co-Variance calculus.
		co- Variance calculus for the equiprobable case : \f$ /frac{1}{N}\sum_{i=0,..,N-1} x_i \cdot y_i - \overline{x}\cdot \overline{y} \f$
*/
		  template <typename T >
		  typename OS_double Covar( const T & x,const T & y) 
		  {
			  return (sum(x)*sum(y))/x.size() - moment1(x)*moment1(y) ;
		  };

/*! 	\brief \f$ \ell^2 \f$ norm.
		 return the square \f$ \ell^2 \f$ norm \f$ \frac{\sum_{i=1,..,N} (a_i)^2 }}{N} \f$.
*/
		  template <typename T >
		  OS_double norm22( const T & Array) 
		  {
			  return sum( Array*Array) ;
		  };
		  template <typename T >
		  bool is_increasing( const T & Array) 
		  {
			  OS_STATIC_CHECK(T::Dim == 1);
			  if ( Array.size() < 2 ) return true;
			  T::const_iterator it_b = Array.begin(),it_e = Array.end();
			  T::value_type local = *it_b;
			  for (++ it_b; it_b != it_e; ++it_b) {
				  if (local >= *it_b) return false;
				  local = *it_b;
			  }
			  return true ;
		  };
/*! 	\brief square \f$ \ell^2 \f$ norm.
		 return the \f$ \ell^2 \f$ norm \f$ \sqrt{\frac{\sum_{i=1,..,N} (a_i)^2 }}{N} } \f$.
*/
		  template <typename T >
		  OS_double norm2( const T & Array) 
		  {
			  return std::sqrt( norm22(Array) );
		  };

		  OS_double moment1( const OS_double &);
//		  template<>
//		  OS_double norm22<particle<OS_double>>( const particle<OS_double> & Array);
		  OS_double sum( const OS_double & );
		  OS_double norm_infty( const OS_double & );


/*! 	\brief square \f$ \ell^1 \f$ norm.
		 return the \f$ \ell^1 \f$ norm \f$ \frac{\sum_{i=1,..,N} |a_i| }}{N}\f$.
*/
		  template <typename T >
		  OS_double norm1( const T & Array)
		  {
			  if (Array.size() == 0) return 0;
			  OS_double local(0.);
			  T::const_iterator it_b = Array.begin(),it_e = Array.end();
			  for (; it_b != it_e; ++it_b) {
				  local += norm1(*it_b);
			  }
			  return local / Array.size();
		  }
		  OS_double norm1( const OS_double & );

/*! 	\brief infimum
		 return the infimum \f$ \inf_i a_i \f$.
*/
		  template <typename T >
		  typename OS_type_helper<T>::value_type inf( const T & Array)
		  {
			  if (Array.size() == 0) return 0;
			  T::value_type result(OS_MAX_REAL);
			  T::const_iterator it_b = Array.begin(),it_e = Array.end();
			  for (; it_b != it_e; ++it_b) {
				  OS_double local_inf = jmmath::inf((*it_b));
				  result = std::min(result,local_inf);
			  }
			  return result ;
		  }
		  template<>
		  OS_double inf<OS_double>( const OS_double & value);

		  template <typename T >
		  OS_double distance( const T & from, const T & to)
		  {
			  return jmmath::norm2(from-to);
		  };
	  }
#endif // jmmath
